SoSe2021
Folienübersicht
Neue Komponenten.
..die es bei mehr als 1 erklärenden Variablen zu berücksichtigen gilt:
- Zusätzliche Datenexploration
- Modellauswahl
- Zusätzliche Modellvalidierung
- Anpassung der Modellgüte
- Interpretation des numerischen R Output und Dummy-Variablen
- Visualisierung der Ergebnisse
Alle X kontinuierlich: Multiple Regression
Multiple Regression
Die Prinzipien der einfachen Regressionsanalyse können auf zwei oder mehr unabhängige Variablen erweitert werden:
Faustregel: nicht mehr als 10 unabhängige Variablen (X)
Voraussetzungen für multiple lineare Regressionsmodelle
- Unabhängigkeit
- Homogenität / homogene Varianz
- Normalität / Normalverteilung
- Linearität
Neu
- Anzahl der Kovariaten < Anzahl der Beobachtungen (p < n)
- Keine lineare Abhängigkeit zwischen den unabhängigen Variablen (keine exakte Multikollinearität)
Problem der Multikollinearität
- Schätzungen von \(\beta\)s und Signifikanztests unzuverlässig.
- Kleine Änderungen in den Daten führen zu extremen Änderungen in den Regressionskoeffizienten.
- Schätzungen zwischen verschiedenen Stichproben der gleichen Grundgesamtheit können sich stark unterscheiden (→ hohe Standardfehler der geschätzten Parameter).
Prüfung
- Paarweise Streudiagramme (mit
pairs()). - Berechnung der Korrelationsmatrix (mit
cor()). - Zusätzliche Bestimmung der Toleranz (TOL) und des VIF-Werts (= ‘Variance Inflation Factor’) für jede Variable.
Variance Inflation Factor
Kurz: VIF
\[VIF = \frac{1}{TOL} = \frac{1}{1-R^2}\]
- \(R^2\) hier aus Modell welches alle X Variablen enthält außer die, für die der VIF-Wert berechnet wird.
- Faustregel: VIF > 3 wird als problematisch eingestuft (manche sagen VIF > 10, hängt auch von der Anzahl der unabhängigen Variablen ab)
- VIF in R:
car::vif()- Wenn X Variablen VIF > 3 haben, Variable mit höchstem VIF entfernen (oder nach biologischem Wissen entscheiden),
- VIFs für verbleibende Variablen neu berechnen,
- und fortfahren, bis alle Variablen VIF < 3 haben.
Variance Inflation Factor
Simulation
Angenommen y ist hauptsächlich abhängig von x1 und nun soll der Steigungskoeffizient \(\beta_1\) geschätzt werden → wenn andere X Variablen, die von x1 abhängen, mit modelliert werden, variiert die Schätzung für \(\beta_1\) bei jeder Wiederholung sehr stark:
Funktionen
# Funktionen
calc_coefs <- function(x1, x2, x3){
y <- x1 + rnorm(length(x1), sd = .3)
# beta für x1 aus 3 Modellen extrahieren
c(coef(lm(y ~ x1))[2],
coef(lm(y ~ x1 + x2))[2],
coef(lm(y ~ x1 + x2 + x3))[2])
}
run_sim <- function(x1, x2, x3,
n_sim = 1000) {
betas <- sapply(1:n_sim,
function(i) calc_coefs(x1, x2, x3))
# Berechnung der Varianz
round(apply(betas, 1, var), 5)
}Vergleich
# Zufällige X Variablen erzeugen n <- 100 x1 <- rnorm(n) x2a <- rnorm(n) x3a <- rnorm(n) # x2 und x3 in Abhängigkeit zu x1: x2b <- x1/sqrt(2) + rnorm(n) /sqrt(2) x3b <- x1*0.95 + rnorm(n)*sqrt(1-0.95^2) run_sim(x1, x2a, x3a)
x1 x1 x1 0.00099 0.00099 0.00099
run_sim(x1, x2b, x3b)
x1 x1 x1 0.00097 0.00229 0.01117
Schätzen und Testen der Koeffizienten
- Auch hier erfolgt die Schätzung über die OLS Methode.
- Y-Achsenabschnitt (Konstante) \(\alpha\):
- Wert von Y wenn alle X Variablen Null betragen → ökologisch weniger relevant
- Partielle Regressionskoeffizienten (Steigungen) \(\beta_p\):
- Einfluss auf Y wenn \(X_p\) sich um eine Einheit erhöht, wobei alle anderen unabhängigen Variablen oder Prädiktoren konstant gehalten werden.
- Separate t-Tests für jeden partiellen Regressionskoeffizienten im Modell: \(H_0: \beta_p = 0\)
Anpassung des Bestimmtheitsmaß
(Multiple) \(R^2\)
\[R^2 = \frac{SS_{Regression}}{SS_{Gesamt}} = 1-\frac{SS_{Residuen}}{SS_{Gesamt}}= 1-\frac{\sum(y_i-\hat{y}_i)^2}{\sum(y_i-\bar{y})^2}\]
- Problem bei multiplen Regression:
- Mit zunehmender Zahl an unabhängigen Variablen steigt die Gesamtstreuung im Nenner, auch wenn diese nicht zur Erklärung beitragen.
- \(\sum(y_i-\bar{y})^2\rightarrow \infty~~~\Rightarrow R^2 \rightarrow 1\)
Adjusted \(R^2\)
\[R_{adj}^2 =1-(1-R^2)\left(\frac{(N-1)}{(N-p-1)}\right)\]
- → d.h. die Freiheitsgrade werden mit berücksichtigt (p = Anzahl der unabhängigen Variablen)
Wie vergleicht und wählt man das beste Modell aus?
- Basierend auf der besten Anpassungsgüte (‘goodness-of-fit’): Wie gut kann ein Modell die Beobachtungen vorhersagen?
- Wie beurteilen wir das?
- Visuell
- Verwendung eines Kriteriums, dass die Diskrepanz zwischen beobachteten Werten und den zu erwartenden Werten des Modells zusammenfasst, z.B.
- erklärte Varianz (\(R^{2}\))
- ‘Akaike’s Information Criterion’ (AIC) und Modifikationen
- ‘Bayesian Information Criterion’ (BIC)
- Ein gutes Buch für diese Thematik wurde von Anderson & Burnham im Jahr 2002 geschrieben: Model Selection and Multimodel Inference - A Practical Information-Theoretic Approach
‘Akaike’s Information Criterion’
AIC()
- Der AIC ist eine relative Messgröße, keine absolute wie \(R^2\).
- Die Werte können nur zwischen Modellen verglichen werden, die an die gleichen Beobachtungen und die gleiche y Variable angepasst wurden.
- →Du kannst keine Modelle vergleichen, die für verschiedene Teilmengen angepasst wurden oder wenn y in Teilen transformiert wurde!
- Die Werte können nur zwischen Modellen verglichen werden, die an die gleichen Beobachtungen und die gleiche y Variable angepasst wurden.
- Je kleiner der AIC, desto besser.
- Faustregel: Wenn der AIC Unterschied zwischen 2 Modellen < 2 ist, dann wähle das einfachere Modell.
- Bei großen Stichproben empfiehlt sich das Bayesian Information Criterion (BIC), da hier der Strafterm ab 8 Variablen ein größeres Gewicht erhält als beim AIC.
Beispiel: airquality
Haben Wind und Sonneneinstrahlung einen Einfluss auf die in New York gemessenen Ozon-Werte?
str(airquality)
'data.frame': 153 obs. of 6 variables: $ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ... $ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ... $ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ... $ Temp : int 67 72 74 62 56 66 65 59 61 69 ... $ Month : int 5 5 5 5 5 5 5 5 5 5 ... $ Day : int 1 2 3 4 5 6 7 8 9 10 ...
# counts of NA purrr::map_int(airquality, ~sum(is.na(.)))
Ozone Solar.R Wind Temp Month Day
37 7 0 0 0 0
air <- airquality %>% tidyr::drop_na(Ozone) %>% tidyr::drop_na(Solar.R)
Beispiel: airquality
Datenexploration
hist(air$Ozone)
air$OzoneS <- sqrt(air$Ozone) hist(air$OzoneS)
pairs(air[ ,c(7,2,3)])
car::vif(lm(OzoneS ~ Wind + Solar.R, data = air))
Wind Solar.R 1.016442 1.016442
Beispiel: airquality
Mögliche Modelle
# Nullmodell: mod0 <- lm(OzoneS ~ 1, data = air) mod1 <- lm(OzoneS ~ Wind, data = air) mod2 <- lm(OzoneS ~ Solar.R, data = air) mod3 <- lm(OzoneS ~ Wind + Solar.R, data = air) # Maximalmodell: mod4 <- lm(OzoneS ~ Wind + Solar.R + Wind:Solar.R, data = air)
Beispiel: airquality
Vergleich des \(R^2\)
summary(mod1)$r.squared summary(mod2)$r.squared summary(mod3)$adj.r.squared # Da 2 X, hier nun das angepasste R^2 auswählen summary(mod4)$adj.r.squared
[1] 0.3703109 [1] 0.162264 [1] 0.4682859 [1] 0.4846256
- Wenn beide erklärenden Variablen im Modell sind, wird deutlich mehr Variabilität erklärt.
- Das Maximalmodell mit der Interaktion bringt aber kaum ein Gewinn an erklärter Variabilität.
Beispiel: airquality
Vergleich des AIC
AIC(mod0, mod1, mod2, mod3, mod4)
df AIC mod0 2 516.0493 mod1 3 466.7085 mod2 3 498.3965 mod3 4 447.8994 mod4 5 445.4023
- Das komplexeste Modell hat den niedrigsten AIC.
- Allerdings beträgt die Differenz zum weniger komplexen Modell nur 2.5 → Wahl fällt auf
mod3.
‘All subsets selection’
- Die eben vorgestellte Methode der Modellauswahl nennt sich auch ‘all subsets selection’.
- Hier werden alle möglichen Teilmengen von Modellen verglichen und das optimale Modell anschl. ausgewählt.
- Nachteile:
- Kann zu einer großen Anzahl von zu vergleichenden Modellen führen.
- Viele Iterationen und zeitaufwendig, wenn p > 5.
- Die Anzahl der Modelle kann jedoch reduziert werden, indem a priori einige wenige zu vergleichende Modelle ausgewählt werden.
- Einige bevorzugen diesen Ansatz, da er einige sehr gute Modelle hervorbringen kann.
Beispiel: airquality
Modellvalidierung 1 - mod3
par(mfrow = c(2,2)) plot(mod3)
Beispiel: airquality
Modellvalidierung 2 - mod3
Zusätzliche Plots: Residuen vs. jede X Variable
air$residuals <- resid(mod3) par(mfrow = c(1,2), mar = c(4,4,1,1)) plot(x = air$Wind, y = air$residuals); abline(h = 0, lty = 2) plot(x = air$Solar.R, y = air$residuals); abline(h = 0, lty = 2)
Beispiel: airquality
Ergebnis
summary(mod3)
Call:
lm(formula = OzoneS ~ Wind + Solar.R, data = air)
Residuals:
Min 1Q Median 3Q Max
-3.9188 -1.3407 -0.1783 1.4120 3.9330
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.241963 0.647547 12.728 < 2e-16 ***
Wind -0.388544 0.048079 -8.081 9.96e-13 ***
Solar.R 0.008855 0.001877 4.719 7.14e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.779 on 108 degrees of freedom
Multiple R-squared: 0.478, Adjusted R-squared: 0.4683
F-statistic: 49.44 on 2 and 108 DF, p-value: 5.706e-16Beispiel: airquality
Visualisierung - Code
Verwende nicht die geom_smooth() Funktion!!
library(modelr)
p_w <- air %>% data_grid(
Wind = seq_range(Wind, 20), # 20 Werte im gleichen Abstand zwischen min und max
Solar.R = mean(Solar.R) # auf Null setzen oder Konstante wie MW wählen
) %>%
add_predictions(mod3) %>%
ggplot(aes(x = Wind)) + ylab("OzoneS") +
geom_line(aes(y = pred), color = "red") +
geom_point(data = air, mapping = aes(y = OzoneS))
p_s <- air %>% data_grid(
Wind = mean(Wind), # diese nun konstant halten
Solar.R = seq_range(Solar.R, 20)
) %>%
add_predictions(mod3) %>%
ggplot(aes(x = Solar.R)) + ylab("OzoneS") +
geom_line(aes(y = pred), color = "blue") +
geom_point(data = air, mapping = aes(y = OzoneS))
gridExtra::grid.arrange(p_w, p_s, ncol = 2)Beispiel: airquality
Visualisierung - Plot
Fragen?
Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License except for the borrowed and mentioned with proper source: statements.
Image on title and end slide: Section of an infrared satellite image showing the Larsen C ice shelf on the Antarctic Peninsula - USGS/NASA Landsat: A Crack of Light in the Polar Dark, Landsat 8 - TIRS, June 17, 2017 (under CC0 license)